#include <stdio.h>
#include <stdlib.h>
#include "mpi.h"
#include <math.h>

double f(double x, double y)
{
  return (-2*(x*(x-1)+y*(y-1)));
}

double max(double x, double y)
{
  return x < y ? y : x;
}

double u_ref(double x, double y)
{
  return x*(x-1)*y*(y-1);
}

int main(int argc, char **argv) {
  if(argc < 2)
    {printf("expected: taille_mat epsilon nb_iteration\n");
      return 0;}
  else
    {
      int max_ite = atoi(argv[3]);
      int N = atoi(argv[1]);
      double epsilon = strtod(argv[2],NULL);
      double h = 1/((double) N-1.0);
      int myrank, tag=99;
      //MPI_Status status;
      MPI_Init( NULL, NULL );
      MPI_Comm_rank( MPI_COMM_WORLD, &myrank );
      if(myrank == 0)
	{
	  int i,j,k;
	  double u[N][N];
	  double unew[N][N];
	  double maximum;
	  //init u
	  for (i=0 ; i<N ;i++)
	    {
	      for (j=0 ; j<N ;j++)
		{
		  u[i][j] = 0;
		}
	    }
	  k=0;
	  do{
	    maximum=0.0;
	    //calcul
	    for (i=1 ; i<N-1 ;i++)
	      {
		for (j=1 ; j<N-1 ;j++)
		  {
		    unew[i][j] = 0.25*(u[i-1][j]+u[i+1][j]+u[i][j-1]+u[i][j+1]-h*h*f((i-1)*h,j*h));
		    maximum=max(maximum,fabs(unew[i][j]-u[i][j]));
		  }
	      }
	    //recopie
	    for (i=1 ; i<N-1 ;i++)
	      {
		for (j=1 ; j<N-1 ;j++)
		  {
		    u[i][j]=unew[i][j];
		  }
	      }
	    k++;
	    //printf("étape %d : epsilon = %f\n",k,maximum);
	  }while(k<max_ite && maximum>=epsilon);
	  int l,m;
	  double diff = 0.0;
	  for (l=1 ; l<N-1 ;l++)
	      {
		for (m=1 ; m<N-1 ;m++)
		  {
		    diff = max(diff,fabs(unew[l][m]-u_ref(l*h,m*h)));
		  }
	      }
	  printf("Différence max avec fonction solution : %f\n",diff);
	}
      MPI_Finalize();
      return 0;
    }
}

